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Abstract 

In this paper, we present a dislocation-density-based three-dimensional continuum model, 
where the dislocation substructures are represented by pairs of dislocation density potential 
functions (DDPFs), denoted by p and p. The slip plane distribution is characterized by 
the contour surfaces of p, while the distribution of dislocation curves on each slip plane 
is identified by the contour curves of p which represents the plastic slip on the slip plane. 
By using DDPFs, we can explicitly write down an evolution equation system, which is 
shown consistent with the underlying discrete dislocation dynamics. The system includes 
i) A constitutive stress rule, which describes how the total stress field is determined in the 
presence of dislocation networks and applied loads; ii) A plastic flow rule, which describes 
how dislocation ensembles evolve. The proposed continuum model is validated through 
comparisons with discrete dislocation dynamics simulation results and experimental data. 
As an application of the proposed model, the “smaller-being-stronger” size effect observed 
in single-crystal micro-pillars is studied. A scaling law for the pillar flow stress cr flow against 
its (non-dimensionalized) size D is derived to be cr flow log (D)/D. 
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1. Introduction 


It is widely agreed that plasticity theories that properly integrate the accumulated knowl¬ 
edge in small-scale physics can facilitate the design of high-end materials. The continuum 
crystal plasticity (CCP) theories (e.g. Rice, 1971; Peirce et ah, 1983; Fleck and Hutchin¬ 
son, 1993; Nix and Gao, 1998; Gurtin, 2002) have shown their values in understanding the 
elasto-plastic behavior of crystals, but they are still phenomenological. On the other hand, 
the (three-dimensional) discrete dislocation dynamical (DDD) models take the dislocation 
microstructural evolution into account based on the fact that plastic deformation of crystals 
is carried out by the motion of a large number of dislocations (e.g. Kubin et al., 1992; Zbib 
et ah, 1998; Fivel et ah, 1998; Ghoniem et ah, 2000; von Blanckenhagen et ah, 2001; Wey- 
gand et ah, 2002; Xiang et ah, 2003; Benzerga et ah, 2004; Quek et ah, 2006; Arsenlis et ah, 
2007; Rao et ah, 2007; El-Awady et ah, 2008; Tang et ah, 2008; Senger et ah, 2008; El-Awady 
et ah, 2009; Chen et ah, 2010; Zhao et ah, 2012; Zhou and LeSar, 2012; Ryu et ah, 2013; Zhu 
et ah, 2013; Zhu and Chapman, 2014b; Wang et ah, 2014). In DDD models, dislocations are 
treated as line singularities embedded into an elastic medium. The kinematics of individual 
dislocations is governed by a collection of laws for dislocation gliding, climb, multiplication, 
annihilation, reaction, etc., and the microstructural changes within crystals are then cap¬ 
tured by the evolution of dislocation curves. DDD models have been well applied to provide 
insights in understanding many plastic deformation processes observed in micro- or nano- 
crystalline structures, such as in thin films and interfaces (e.g. von Blanckenhagen et ah, 
2001; Weygand et ah, 2002; Quek et ah, 2006; Zhou and LeSar, 2012; Wang et ah, 2014) and 
in micro-pillars (e.g. Rao et ah, 2007; El-Awady et ah, 2008; Tang et ah, 2008; Senger et ah, 
2008; El-Awady et ah, 2009; Ryu et ah, 2013). However, three-dimensional DDD models 
become too computationally intensive when the specimen size exceeds the order of several 
microns. 

Therefore, a successful dislocation-density-based theory of plasticity (DDBTP) whose 
associated length scale lies between CCP’s and DDD’s is still highly expected. The devel¬ 
opment of DDBTP dates back to the works of Nye (1953), where a dislocation network is 
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represented by a continuously distributed second-order tensor, known as the Nye disloca¬ 
tion tensor. Nowadays with more knowledge in physics taking place on smaller scales, a 
successful DDBTP should be constituted by laws that are consistent with the underlying 
discrete dislocation dynamics from the following two aspects: i) A constitutive stress rule 
to determine the stress field in the presence of a continuous dislocation density distribution 
and applied loads; ii) A plastic flow rule to capture the motion of dislocation ensembles (in 
response to the calculated stress field), which results in plastic flows in crystals. 

As the simplest dislocation configuration, systems of straight and mutually parallel dis¬ 
locations have been analyzed relatively well at the continuum level (e.g. Groma et al., 2003; 
Berdichevsky, 2006; Voskoboinikov et ah, 2007; Kochmann and Le, 2008; Hall, 2011; Oz- 
top et ah, 2013; Geers et ah, 2013; Le and Guenther, 2014; Schulz et ah, 2014; Zhu and 
Chapman, 2014a; Le and Guenther, 2015). In this case, each dislocation can be treated as 
a point singularity in a plane that is perpendicular to all dislocations. As a result, the Nye 
dislocation density tensor is reduced to several scalar dislocation density functions. The 
geometric complexity of the dislocation networks is dramatically reduced in this case. How¬ 
ever, the development of three-dimensional DDBTP is still far from satisfactory despite a 
number of valuable works (e.g. Nye, 1953; Kroener, 1963; Kosevich, 1979; Nelson and Toner, 
1981; Mura, 1987; Head et ah, 1993; El-Azab, 2000; Acharya, 2001; Svendsen, 2002; Arsenlis 
and Parks, 2002; Sedlacek et ah, 2003; Alankar et ah, 2011; Sandfeld et ah, 2011; Engels 
et ah, 2012; Hochrainer et ah, 2014; Li et ah, 2014; Cheng et ah, 2014). One of the main 
barriers in establishing a successful three-dimensional theory is due to the fact that the com¬ 
plex networks of curved dislocation substructures make the upscaling of discrete dislocation 
dynamics extremely difficult. 

To overcome such difficulties, Xiang (2009) introduced the idea of a coarse-grained dis- 
registry function (CGDF), which is defined to approximate the exact disregistry function 
(plastic slip) used in the Peierls-Nabarro models (Peierls, 1940; Nabarro, 1947; Xiang et ah, 
2008), by a smoothly varying profile without resolving details of dislocation cores. By this 
way, the density distribution of a discrete curved dislocation network in a single slip plane 
(after local homogenization) can be simply represented by the scalar CGDF (more precisely, 
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the spatial derivatives of CGDF), and dislocation dynamics on the slip plane at the contin¬ 
uum level is explicitly formulated in terms of the evolution of the CGDF (Zhu and Xiang, 
2010). Using this representation of CGDF for dislocation density distribution, connectivity 
condition of dislocations is automatically satisfied. The underlying topological changes of 
dislocations are automatically handled by the evolution equation of the CGDF, and no law 
for dislocation annihilation needs to be further imposed. The dynamics of dislocations in the 
continuum model is derived from the DDD model, and the dislocation velocity in the contin¬ 
uum model depends on a continuum version of the Peach-Koehler force on dislocations. It 
has been rigorously shown by Xiang (2009) that in the continuum model, the Peach-Koehler 
force due to the resolved shear stress of a family of curved dislocations can be decomposed 
into a long-range dislocation interaction force and a short-range self line tension force, and 
they can both be expressed in terms of the spatial derivatives of CGDFs. The Frank-Read 
source, which is one of the major mechanisms for dislocation multiplication, is also well 
incorporated into this continuum framework (Zhu et al., 2014). As one application of this 
continuum model using CGDFs, a two-dimensional Hall-Petch law, which relates the flow 
stress of a polycrystal not only to the physical dimension of its constituent grains, but also 
to the grain aspect ratio, is derived without any adjustable parameters (Zhu et ah, 2014). 

In this paper, we generalize our previous single-slip-plane model to that for dislocation 
ensembles in three-dimensions, where the density distribution of dislocations is locally co¬ 
determined by an in-plane dislocation density distribution and a slip plane distribution. To 
take into account the spatial variation from these two aspects, we define a pair of dislo¬ 
cation density potential functions (DDPFs) for each active slip system. One DDPF 
0 is employed to describe the slip plane distribution (after local homogenization) by its 
contour surfaces, and the other DDPF 0 is defined such that 0 restricted on each slip plane 
describes the plastic slip across the slip plane and identifies the density distribution of dislo¬ 
cation curves (after local homogenization) on that plane. Here we name 0 and 0 by density 
potential functions, because the Nye dislocation density tensor is represented in terms of 
the spatial derivatives of these two functions. As our previous continuum model in a single 

slip plane (Xiang, 2009), the major advantage of this three-dimensional continuum model 
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lies in its simple representation of distributions of curved dislocations (after local homoge¬ 
nization) using two scalar DDPFs, which automatically satisfies the connectivity condition 
of dislocations. 

To derive the constitutive stress rule in the continuum framework with the DDPFs, 
we sequentially express the Nye dislocation density tensor, the plastic distortion and the 
elastic strain tensor in terms of the DDPFs. As in our previous continuum model in a 
single slip plane (Xiang, 2009), the continuum Peach-Koehler force due to the resolved 
shear stress consists of a long-range dislocation interaction force and a short-range self line 
tension force. The long-range stress field is determined by the derived constitutive stress 
rule and the equilibrium equations along with boundary conditions, and a finite element 
(FE) formulation is proposed to compute this long-range stress field. The local self line 
tension effect can be explicitly formulated in terms of the spatial derivatives of 0 and -0. 
The plastic flow rule is described by evolution equations of the DDPFs 0 and 0. For face- 
centered-cubic (fee) crystals considered in this paper, the motion of dislocations is limited 
to their respective slip planes, and the plastic flow rule is given by an evolution equation 
of 0. Frank-Read sources in our single-slip-plane continuum model (Zhu et al., 2014) is 
generalized to three-dimensional case with continuous distributions of a number of sources. 
The derived equations form a closed system evolving in time as summarized in Eqs. (54) to 
(57) in Sec. 2.9. 

With this continuum model, we investigate the size effect on crystalline strength widely 
observed in the uniaxial compression tests of monocrystalline micro-pillars (e.g. Uchic et al., 
2004, 2009; Jang et al., 2012). Practically, an empirical power law is adopted to relate 
the pillar flow stress <7fl ow to the pillar size D by crfl ow ~ D~ m , where m is found to be 
from 0 to 1, varying by study (Uchic et al., 2009). Typically there are two classes of 
models proposed to rationalize this size effect. The first type falls into the family of the 
“dislocation starvation’’ models (Greer et al., 2005; Greer and Nix, 2006). They argued 
that a crystal small in size does not provide enough space for dislocation multiplication 
and the flow strength gets increased as a result. The second category of models attribute 

the observed size effect to the stochastics of dislocation source lengths in the small-size 
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specimens (Parthasarathy et al., 2007). There are also models using statistical approaches 
to reproduce the power law expression (e.g. Gu and Ngan, 2013). Analysis of heterogeneous 
deformation in single crystal micropillars under compression has been performed by using 


a hybrid elasto-viscoplastic simulation model which couples DDD model, and a scaling law 
(jfjow ~ .D~"log.D + aD~ l with two parameters a and 0 < n < 1 was proposed (Akarapu 


et al., 2010). In the last part of this paper, we apply the proposed continuum model to 
study the size effect on the strength of micro-pillars, by following the trace of the stochastic 
source length models which have been employed to rationalize the size effect observed in 
DDD simulations (e.g. El-Awady et al., 2008; Senger et ah, 2008; El-Awady et ah, 2009; 
Shao et ah, 2014). We find that the flow stress scales with the sample size by 



( 1 ) 


This relation is validated through comparison with experimental data conducted in several 
fee crystals for pillar size ranging from submicrons to tens of microns. 

The rest of this paper is organized as follows. In Sec. 2, the (three-dimensional) DDPFs 
are introduced, and this is followed by the derivation of the constitutive stress rule and the 
plastic flow rule needed at the continuum level. In Sec. 3, numerical schemes for the derived 
equation system are presented. In Sec. 4, the derived continuum model is validated through 
comparisons with DDD simulation results. In Sec. 5, the continuum model is applied to study 
the size effect arising in the uniaxial compression tests of single crystalline micro-pillars. 

To better illustrate the derivation of the continuum model, following notations are used 
throughout the article unless specified. The Cartesian coordinates are denoted by r = 
(x,y,z). The i-th entry of a vector, for example r, is denoted by ?y, and the ij-th entry 
of a second-order tensor, for example cr, is denoted by Unless specified, the following 
notations are used: the vector gradient (Vu)^ = dui/drj ; the cross product (m x n)* = 
Yj k=i e ijk m jnk with Eijk the permutation tensor; the inner product of two vectors m • n = 
Y^=i m i n u ^ ie i nner product of two second-order tensors a : (3 = Yl 3 =\ a ijPij'i the inner 
product of a fourth-order tensor and a second-order tensor C : (3 — Yt i=i the 

magnitude of a vector |u| = • u; the symmetric part of a second-order tensor sym(o:) = 


6 



(a + a T )/ 2; the outer product of two vectors (a 0 b)^- = atbj ; and the row “curl” of a 
second-order tensor (V x ct)^ = e jki®u,k- 

2. Continuum plasticity model based on dislocation density potential functions 

In this section, we first review the continuum model for dislocation dynamics in one slip 
plane using a two-dimensional coarse-grained disregistry function (Xiang, 2009; Zhu and 
Xiang, 2010; Zhu et ah, 2014). Then by using the DDPFs, we build the three-dimensional 
continuum model for dislocation dynamics and plasticity. 

2.1. Review of the continuum model in a single slip plane 

We consider describing a given discrete dislocation network in a single slip plane and 
with the same Burgers vector b (e.g. the configuration on the top left of Fig. 1(a)) by 
a dislocation continuum. For this purpose, we use two held quantities at the continuum 



Figure 1: (a) Turning a discrete dislocation network in a single slip plane into a dislocation continuum. 
The discrete dislocation network with Burgers vector b is approximated by a family of smoothly varying 
dislocation curves with local line direction l and dislocation spacing c?i n , averaged locally over a representative 
area IF of the discrete dislocation network, see Eq. (2). (b) Representation of the dislocation continuum in 
(a) (on the top right) using the CGDF ^ 2 d- The i-th dislocation curve corresponds to the contour of b 2 d 
with height </> 2 d = ib. 
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level: the average dislocation line direction l and the average dislocation spacing d m , which 
equals the reciprocal of the net dislocation length per area /?g d . The values of the two held 
quantities l and d m at each point in the continuum model come from the averaging over a 
representative area IT of size e in the discrete model, where d- m « e « D with D being 
the sample size in the continuum model. Under this assumption, all dislocations inside fT 
can be treated as line segments, as schematically shown in Fig. 1(a). By superpositioning 
all dislocation segments inside Q e , we can obtain a super line segment denoted by L, i.e. 
L = s i f° r line segment s,; e O e . The quantities of interest at the continuum level can 
then be defined by 


l = lim 


L 

1 

= lim 

|1T| 



L 

P 2d 

e/D—^O 

|L| 


( 2 ) 


where |fF| is the area of IT. 

By performing such averaging process everywhere, the original discrete dislocation net¬ 
work is turned into a dislocation continuum as schematically shown in Fig. 1(a) (on the 
top right) and we can then use held variables to express it. Xiang (2009) introduced a 
coarse-grained disregistry function (f> 2 d, such that the ?'-th dislocation curve in the averaged 
dislocation continuum corresponds to the contour of 0 2 d with height 0 2 d = ib, see Fig. 1(b). 
With the function 0 2 d, the local dislocation line direction and inter-dislocation spacing can 
be given by 


and 


. _ 1 f d(j) 2d dcj)2d \ 

|V 2d 02d| V dy dx ) 


/ _ 1 _ b 

m P 2 g d |V 2d 0 2d r 


( 3 ) 


( 4 ) 


respectively, where V 2d = ( ^ 


A fact used to derive these formulas is that the normal 
vector n of the i-th dislocation curve in the averaged dislocation continuum 


V2d02d 
| V2d02d| ’ 


( 5 ) 


which is the unit vector in the slip plane and normal to l. 








With the introduction of 0 2 d, we can capture the resolved shear stress including that due 


to the long-range dislocation interaction by an integral 


2d _ _h_ 
long 4 t r 


x 


x 


94>2d(x,y) 

dx 


+ (y-y) 


9<t> 2a(x,y) 

9y 


((x — x) 2 + (y — |/) 2 ) 3 / 2 


-dxdy 


/iv 

^ 47t(1 — u)b 2 


( ^ 9<f>2d{x,y) 
l 1 dx 


+ b 2 d(p2d d { p y) ^j {bi(x -x) + b 2 (y - y)) 
((x — x) 2 + (y — y ) 2 ) 3 / 2 


dxdi/, 


( 6 ) 


where fj, and v are the shear modulus and the Poisson’s ratio, respectively, and a contribution 
due to the local line tension effect 


2d _ y bK 
Tself 4 t r 


( 

1 + 19 
l-V 


\ 


5 fe‘ + ^) 2 /6 2 |, 

- log 


' f)rhr> j ' 



\ 

+ 1 

/ 


(7) 


where r c is a parameter depending on the dislocation core and k = -V - (V0 2 d/| V0 2 d|) is the 
local (average) curvature of the dislocation. These stress formulas in the continuum model 
were derived rigorously from the discrete dislocation model by asymptotic analysis (Xiang, 


2009). 


In the single-slip-plane continuum model, 0 2( j measures the plastic slip across the slip 


plane in the direction of the Burgers vector, and the plastic flow is governed by an evolution 
equation of 0 2 a (Zhu and Xiang, 2010; Zhu et ah, 2014): 


d02d 

dt 



+ 



S2d- 


( 8 ) 


Here v n is the dislocation moving speed along its normal direction (Zhu and Xiang, 2010) 


v„ = m,(rS g + T™,)b, 


(9) 


where m g is the dislocation glide mobility and the applied stress can also be included in 
the above equation. The right-hand term s 2( j in Eq. (8) formulates the effect due to the 
dislocation multiplication by Frank-Read sources (Zhu et al., 2014), which will be reviewed 
in detail in Sec. 2.8.1. 
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2.2. Dislocation substructures in three dimensions represented by dislocation density poten¬ 
tial functions 

From now on, we present our continuum model for the dynamics of dislocation struc¬ 
tures in three dimensions. In this subsection, we introduce the representation of dislocation 
substructures in three dimensions using DDPFs. 




(a) (b) 

Figure 2: (a) A discrete dislocation network in three dimensions is approximated by a dislocation continuum, 
after local homogenization of dislocation ensembles within some representative volume Q e . These dislocations 
in the network are associated with the same slip system with Burgers vector b and slip plane normal direction 
m. The dislocation continuum is characterized by the dislocation line direction Z, the slip plane spacing d s \ 
and the in-plane dislocation spacing di n averaged over the representative volume Q e , see Eqs. (11) and (12). 
(b) Representation of the dislocation continuum in (a) using a pair of DDPFs </> and ip. The DDPF ip is 
employed such that the j -th slip plane is the contour plane of ip of height ip = jb. The DDPF <p is defined 
such that <p restricted on each slip plane describes the dislocation distribution on that plane, i.e., the i-th 
dislocation on the slip plane is the contour line <p = ib. 

We first focus on a given discrete dislocation network in a single slip system with Burgers 
vector b and slip plane normal direction m. We take a representative cuboid volume with 
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size e, the Nye dislocation density tensor over D e is (Nye, 1953; Kroener, 1963; Kosevich, 
1979) ot = JTb® Sj/|h2 e | for dislocation line segment s* G D e , where 10 e | is the volume of 
Q e , see Fig. 2(a). Denoting L = s ?: for all dislocation line segments Sj e O e which is the 

superposition of all dislocation segments inside fP, and assuming that e « D where D is 
the domain size of the continuum model, we write the Nye dislocation density tensor over 
IP as 

ol = p g b <8) Z, (10) 


where Z is the average dislocation line direction over O e and p g is the dislocation density 


over fZ e (net length per unit volume): 


Z 


lim 777, 

e/D-sU L 


lim — 

e/D ->0 |fZ e 


( 11 ) 


See Fig. 2(a) for an illustration of this average. We characterize p g by two variables 


Pg 


dsl d] n 


( 12 ) 


where d s i is the slip plane spacing and d in is the in-plane dislocation spacing, see Fig. 2(a). 

By performing such averaging process everywhere, the original three-dimensional discrete 
dislocation network is turned into a dislocation continuum in three dimensions as schemati¬ 
cally shown in Fig. 2(a). We introduce a pair of dislocation density potential functions 
(DDPFs) (j) and to represent this resulting dislocation continuum consisting of dislocations 
in the same slip system. The DDPF -0 is employed to describe the distribution of the slip 
planes in the averaged dislocation continuum, such that the j'-th slip plane is the contour 
plane of 0 with 0 = jb. In this paper, we focus on fee crystals, in which dislocations stay in 
their respective slip planes. We further assume a uniform distribution of the slip planes with 
prescribed slip plane spacing cZ s i- Under these conditions, the DDPF 0 for the slip system 
with Burgers vector b and slip normal direction m can be written as 


0(r) 


6m • (r — r°) 

da 


( 13 ) 


where r° is some point on the 0-th slip plane, see Fig. 2(b). 
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With slip planes in the averaged dislocation continuum described by 0, we introduce 
another DDPF 0 to describe the distribution of dislocations on each slip plane, such that 
0 restricted on each slip plane is the two-dimensional CGDF 0 2 d in the single-slip plane 
model reviewed in the previous subsection. That is, the i-th dislocation curve on the j'-th 
slip plane in the averaged dislocation continuum can be represented by 


(r|0(r) = ib, 0(r) = jb , for integers i,j}, (14) 

see Fig. 2(b). As 0 2 d in the single-slip plane model, here 0 restricted on each slip plane 
describes the plastic slip across the slip plane in the direction of the Burgers vector. Local 
geometric quantities of dislocations and the Nye dislocation density tensor in the averaged 
dislocation continuum are simply expressed in terms of these DDPFs 0 and 0, see the 
following subsections. 

The major advantage of this three-dimensional continuum model characterized by DDPFs 
lies in its simple representation of distributions of curved dislocations (after local homog¬ 
enization). This representation also automatically satisfies the connectivity condition of 
dislocations, see Sec. 2.4. In addition, dislocation annihilation within the same slip plane is 
automatically handled as in the previous single-slip-plane model (Xiang, 2009; Zhu and Xi¬ 
ang, 2010). When there are multiple slip systems activated, one can assign a pair of DDPFs 
to each active slip system. 

In body-centered-cubic (bcc) crystals or fee crystals at high temperatures where dislo¬ 
cations are able to move out of their slip planes by climb, dislocations no longer stay on 
planar slip planes. Dislocation networks in these crystals can still be represented under the 
continuum framework characterized by the DDPFs 0 and 0. In these cases, the contour 
surface 0(r) = jb may describe a curved surface rather than a plane, and the local slip plane 
normal direction in the averaged dislocation continuum is determined by 

V0 


m = 


IV-01 


(15) 


2.3. Geometrical structures of dislocation continuum described by DDPFs 


Now we present expressions of the local geometric quantities of dislocations in the aver¬ 
aged dislocation continuum in terms of the DDPFs 0 and 0. The slip plane normal direction 
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is given by Eq. (15). The local dislocation line direction l is 


_ V0 x V-0 
| V0 x V-01 


( 16 ) 


This is because the dislocation is contained in both the contour surfaces of 0 — ib and if = jb 
(see Eq. (14)), thus it is perpendicular to both V0 and V0. The local dislocation normal n 
(which is in the slip plane) is also calculated in terms of DDPFs by 


V0 x (V0 x V0) 
|V0 x V-011V-01 


(17) 


Moreover, by the Frenet-Serret formulas, we have kw = (/ ■ V/i, / ■ V/ 2 , l • V/ 3 ), where ft is 
the curvature of the local dislocation. Thus with the expressions for l and n in Eqs. (16) 
and (17), k is also represented by the DDPFs as 


3 

k = n • (fen) = ■ VZj. (18) 

i= 1 

Finally, the in-plane dislocation spacing and the slip plane spacing are respectively 


d it 1 


| V0 x m| ’ 


del = 


IV# 


(19) 


Here we have used the definition in Eq. (14) and the fact that the gradient of 0 over the slip 
plane is V in _ p i ane 0 = V0 — (m • V0)m = V0 x m. Note that the above formulas hold for a 
general DDPF 0 for the definition of dislocations in Eq. (14), including but not limited to 
the particular expression in Eq. (13). 


2-4- Representation of the Nye dislocation density tensor by DDPFs 

After the local homogenization process described in Sec. 2.2, the Nye dislocation den¬ 
sity tensor at. for a single slip system is calculated by Eq. (10) in terms of the averaged 
dislocation line direction / and the net length per volume dislocation density p g over some 
representative volume, and p g is given by Eq. (12) using the slip plane spacing d s i and the 
in-plane dislocation spacing d m . Further using the formulas in Eqs. (15), (16) and (19) with 
the DDPFs 0 and if, we have 


Pg 


|V0 x m| ■ |V0| 


|V0 x V0| 


b 2 


13 


h 2 


( 20 ) 








Therefore, incorporating Eqs. (16) and (20) into Eq. (10) gives the following expression 
for the Nye dislocation density tensor based on the DDPFs 


a = — <g> (V0 x Vxji) 
b z 


/ / ( chp_ dcj)_ chp_ d(f) 

1 \ f)ni r)n r)y 


1 

62 


dz dy dy dz 

h ( _ dp dp 

2 ydz dy dy dz 

\b. 


dp dp 
dz dy 


dtl> d<j> 
dy dz 


u ( dtp 0p a dt/) d<j> \ 

1 v Ox dz dz dx ) 

h\ 

( dip dip 

dp dp 

p dy dx 

dx dy 

i / di/j dp dip dip \ 

2 V dx dz dz Ox ) 

b 2 1 

( dip dip 

dp dp 

p dy dx 

Ox dy 

i f dip dip dp) dip \ 

3 V Ox dz dz Ox ) 

63 1 

( dp dp 

dp dp 

l dy dx 

Ox dy 


( 21 ) 


For the case with multiple slip systems, we use one pair of DDPFs 0 A and ^ A for the 
A-th slip system, and the Nye dislocation density tensor is expressed by 

b A 


ct = 


E 


(fr 


A)2 


<g> (V0 A x V*/> A ). 


It is easy to check that a. given by Eq. (22) satisfies 

V ■ OL = 0, 


( 22 ) 


(23) 


which means = 0 for a H * = 1,2 and 3. This is the connectivity condition of 

dislocations, which is due to the fact that dislocation lines can never begin or end inside the 
sample. 


2.5. Constitutive stress rule 

In this subsection, we present the constitutive stress rule in our continuum plasticity 
model using DDPFs. The derivation is based on the constitutive relations commonly used in 
classical dislocation-density-based models, including (when the specimen experiences small 
deformations): 

• The total distortion, which is the gradient of the displacement field u, can be decom¬ 
posed into an elastic distortion (3 C and a plastic distortion (3 P 

Vu = /3 C + /3 P . (24) 


The Nye dislocation density tensor equals the curl gradient of the plastic distortion 


Vx/3 p = -a = -^ 


(b 


\\2 


(V0 A x V^ A ). 


(25) 
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( 26 ) 


• The stress field cr satisfies the Hooke’s law (e.g. the isotropic case): 

<r = 2/ie e + y^-tr(e e )I. 

where e e is the elastic strain tensor related to the elastic distortion by e e = sym (/3 e ), 
tr(e e ) is the trace of e e , and I is the 3x3 identity matrix. 


The equilibrium condition in absence of body forces is 


V • cr = 0. 


(27) 


By generalizing the expression for f3 p for the case of discrete dislocations (e.g., Sec. 1.6 
of Mura, 1987), and using the fact that </> describes the plastic slip across the slip plane in 
the direction of the Burgers vector and if) describes the distribution of the slip planes, in the 
continuum model we have 


/3 p = -£ 


(fr 


M2 


(b A (g) v^ A ). 


(28) 


It can be checked that /3 P given by Eq. (28) satisfies Eq. (25). A comparison between Eq. (28) 
and its counterpart in continuum crystal plasticity theories suggests that the scalar 

cj) x |V^ A | 


7 A = 


b x 


(29) 


measures the magnitude of the plastic shearing in the direction of Burgers vector b A in the 
A-th slip system. 

By using Eqs. (24) and (28), we can express the elastic strain tensor by 


e e = sym(Vu) + ^ ^ A ^ 2 sym(b A <g) V0 A ) 


(30) 


Incorporating Eq. (30) with the Hooke’s law (26) and using the fact that b A • V?/> A = 0 (the 
Burgers vector is always perpendicular to the normal direction of the slip plane), we obtain 


cr = £ : Vu + 2/i ^ ^sym(b A <g> V^ A ), 
where the symmetric fourth-order tensor C is defined such that 


(V-u)I 


(31) 


C : Vu = 2/i ( sym (Vu) + 

15 


1 - 2u 


(32) 








When a solid body is purely elastic (dislocation free), (f) x = 0 and Eq. (31) becomes 
<t = C : Vu, which is exactly the constitutive relation in the linear elasticity theory. 

In summary, with given dislocation distributions described by DDPF pairs <f x and if> x 
for A varies over all the slip systems, the stress held er is obtained by solving Eq. (31) and 
the equilibrium condition in Eq. (27) over the sample ff, subject to appropriate boundary 
conditions on the boundary <9Q. One boundary condition is the displacement boundary 
condition imposed on part of the boundary denoted by <9h2d ; 

u|dn d = u b , (33) 

another boundary condition is the traction boundary conditions imposed on the rest of the 
boundary denoted by <9f2 t : 

o-|an t • k = t b , (34) 

with k the outer unit normal to the surface dfl t - Here dfl = U <9A2f 


2.6. Short-range dislocation line tension effect 

It is well-known that there are various short-range dislocation interactions that play 
important roles in the plastic deformation processes, in addition to the long-range effect 
of dislocations. In this subsection, we first show that the continuum stress formulation 
presented in the previous subsection is for the long-range dislocation interaction. We then 
present the formulation for the short-range line tension effect in the DDPFs based three 
dimensional continuum model. 

When the medium is the whole three dimensional space M 3 , using the Nye dislocation 
density tensor formula in Eq. (22) in our continuum model, the stress formulas given by the 
DDD model and the continuum model are as follows (Hirth and Lothe, 1982; Mura, 1987). 
In the DDD model, the stress held is 


cr dd( r ) = 


2vr 


n x , 

A i =1 

_h 

47t(1 


sym 


n x 


b A x (r 


— r|3 


Ids 


bjEE/, (Z ■ (b A x V))(V ® V — IV' 


( 35 ) 


)|r — ijds, 
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where 7 * denotes the i-th dislocation in the A-th slip system, N\ is the number of dislocations 
of the A-th slip system, r goes over all points on 7 *, s is the arclength of 7 A, and V denotes 
the gradient with respect to r. In the continuum framework characterized by the DDPFs, 
the stress field is 

CT ( r ) = E S - Vln ( b ® ( V ^( f ) x V^ x (f))i AV 

A . (36) 

+ E 4t(6 a) 2 (1 _ v) Jj bX ■ V^ A (f))(V^(f) ■ V)(V 0 V - IV 2 )|r - f |dV, 

where dV is an infinitesimal volume associated with f. Since the Nye dislocation density of 
the DDD model is averaged into the Nye dislocation density of the continuum model in the 
coarse-graining process as described in Secs. 2.2 and 2.4, from Eqs. (35) and (36), we have 

Odd —» rr (37) 

in the coarse-graining process from the DDD model to the continuum model. 

When we consider a material with finite size, an image stress is added in the DDD 
model to accommodate the boundary conditions (Van der Giessen and Needleman, 1995). 
Following the same coarse-graining process, such a DDD stress solution gives a stress solution 
that satisfies all the equations and boundary conditions of the continuum model, which is 
the same as that obtained by directly solving the continuum model due to the uniqueness 
of the solution. Therefore in this case, the stress held solved from the continuum model is 
also the leading order approximation of the stress held of the DDD model. 

The stress held cr solved from the elasticity system in the continuum model describes the 
long-range elastic interaction of dislocations. There are also various short-range dislocation 
interactions that play important roles in the plastic processes. How to capture these short- 
range effects at the continuum level is an important issue in the development of dislocation 
based plasticity theories. 

One way to incorporate the short-range dislocation interactions in the continuum model 
is to add complementary resolved shear stresses that describes these effects, because the 
dynamics of dislocations depends on the resolved component of the stress held in the slip 
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plane (Hirth and Lothe, 1982). For the A-th slip system, the resolved shear stress due to 
the long-range dislocation interaction in the continuum model is 

V0 A 


r A 

long 


|V0 ; 


(38) 


Recall that m A = V0 A /|V0 A | is the normal direction of the slip plane. The principle to 
include additional shear stresses in the continuum model is to obtain a better approximation 
to the shear stress in the DDD model 


A _> _ A _ __ A | _ A 

' dd ' 'total 'long ' '*= 


' total 


short 


(39) 


in the coarse-graining process from the DDD model to the continuum model, where r s ) lort 
includes contributions to the resolved shear stress due to all the important short-range 
effects, and r dd = ^ • cr <m • m A is the resolved shear stress of the A-th slip system using the 
DDD model. 

In this paper, we focus on the dislocation line tension effect which is one of the important 
short-range dislocation effects, and the total shear stress in the continuum model is 


Rotai = hong + Rein (40) 

where r s A lf is the contribution to the resolved stress due to the line tension effect. The dis¬ 
location line tension effect plays crucial roles for example in dislocation multiplication by 
Frank-Read sources. For dislocation distributions in a single slip plane, we have used asymp¬ 
totic analysis to rigorously derive the expression of the line tension effect in the continuum 
model from the DDD model (Xiang, 2009) as reviewed in Sec. 2.1. Here for dislocation 
distributions in three dimensions represented by DDPFs, for the A-th slip system, r s A lf is 
given by 

_ fib x fl + v 3n |V0 A | 2 (b A ■ V0 A ) 2 \ / 6 A |V0 A | 

Tself “ 4 tT \1 - v ~ 1 - ‘ (6 A ) 2 |V0 A x V0 A | 2 ) K ' ° S V2vrr c |V^ A x V0 A 

where the curvature k of the local dislocation is calculated by Eq. (18). In order to derive 

this formula, we have used the fact that the DDPF 0 restricted onto each slip plane which 

is a contour surface of the DDPF -0 reduces to our previous single slip plane model as well 

as the expressions of geometric characterizations of the local dislocation given in Sec. 2.3. 
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There are other short-range effects to be included in r s hort, such as that due to the short- 
range interaction between dislocations from different slip systems. These will be discussed 
elsewhere. 


2. 7. Plastic flow rule 

In our continuum model characterized by DDPFs, the plastic flow rule is given by 


■\ a |V0 a xV^ 


|V^ 


= s 


(42) 


for the A-th slip system, where <p x = d(p x /dt, v x is the speed of the local dislocation in its 
normal direction in the slip plane, and s a formulates the dislocation generation by Frank- 
Read sources to be discussed in details in the next subsection. Being the three-dimensional 
version of Eq. (8), Eq. (42) is also established based on the conservation of plastic shear 
slips (Zhu and Xiang (2010), see also the level set DDD method in Xiang et ah (2003)) and 
the fact that the gradient of fl x over the slip plane is V</> A x m A , where m A = V r fl x /\ V^ A | 
is the normal direction of the slip planes. Note that there is no need to assign extra rule 
for dislocation annihilation on the same slip plane, which is automatically handled by the 
topological changes in the contours of the DDPF 0 A during its evolution. 

The dislocation speed u A in Eq. (42) is determined by a mobility law following the DDD 
models: 

u A = m g b x r x tal = m g b x (r? ong + r A lf ) , (43) 

where m g is the dislocation glide mobility, is the component of the long-range stress 
field er determined by the elasticity problem in Sec. 2.5 and resolved in the A-th slip system 
given in Eq. (38), and r A lf is the effective stress due to the line tension effect given in Eq. (41) 
in the previous subsection. 

In general, boundary conditions are needed for Eq. (42). One extreme case is that 
dislocations can exit the specimen freely, which can be described by the Neumann boundary 
condition = 0 where k s is the outer normal direction of the intersection of the slip 
plane and the specimen surface (k s = m x (m x k)/|m x (m x k)| where k is the outer 


19 



normal direction of the specimen surface). Another extreme case is that the dislocations are 
impenetrable to a specimen surface, where <f x is fixed on the specimen surface. 


The rate of Nye dislocation density tensor a, the rate of the plastic distortion (3 P and 
the rate of the scalar plastic shear j x can all be easily calculated using <p x in Eq. (42) and 
the expressions of ck, /3 p and y A in terms of </> A and ip x in Eqs. (22), (28) and (29). Note 
that in this paper, the DDPF if x that describes the distribution of the slip planes is fixed. 


Many other quantities that are useful in understanding the plastic behavior of crystals 
can also be expressed in terms of the DDPFs </> and ip. For example, the total dislocation 
density (length per unit volume) within the specimen 12 can be formulated by 



(44) 


where p 0 is the dislocation density due to the initial distribution of dislocation segments of 
the Frank-Read sources (see the next subsection). Moreover, the total plastic strain rate, 
which is conventionally defined to be the rate of area swept by all dislocations multiplied by 
the length of the respective Burgers vector per volume, can be expressed by 



(45) 


2.8. Incorporation of Frank-Read sources into the continuum model 

In this subsection, we present the expression for the source term s x in the plastic flow 


rule in Eq. (42) for a continuous distribution of Frank-Read sources in three dimensions, 
which is generalized from our previous continuum model for individual sources on a single 
slip plane. 

2.8.1. Review of Frank-Read sources in the single-slip plane continuum model 

We first review our continuum model for individual sources on a single slip plane (Zhu 
et ah, 2014). 

A Frank-Read source is a dislocation segment pinned at its two ends. When the resolved 
shear stress r acting on it exceeds a critical value, known as the activation stress r c , it will 
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keep injecting dislocation loops to the system (Hirth and Lothe, 1982). The time it takes 
for a Frank-Read source to perform an operating cycle is known as the nucleation time t nuc . 

However, if observed at the continuum level, what we see is not the detailed loop-releasing 
process, but continuous dislocation flux originating from a small source region. In the single¬ 
slip plane model (Zhu et ah, 2014) reviewed in Sec. 2.1, the operation of a Frank-Read source 
is controlled by three parameters all determined by the underlying DDD model: the source 
activation stress r c , the source operating rate which equals to l/f nuc , and the source region 
h2g d which is the two-dimensional region enclosed by a newly released dislocation loop. 

The nucleation stress r c is evaluated by adopting the critical stress formula given by 

Tc = iff log (l) • (46) 

where C s depends on the source character and the Poisson’s ratio v (with v — 1/3, C s — 1 
for an edge-oriented source and C s = 1.5 for an screw-oriented source), l is the length of the 
Frank-Read source, and r c is a parameter depending on the dislocation core. 

The nucleation time t nuc is calculated to be 


t 


nuc 


Qchl 

m gK\ r \ - TcV 


(47) 


where r is the resolved shear stress due to the long-range stress field, and Q c h depends 
only on the source orientation fitted from the DDD simulation (Q c h = 6.1278 for edge- 
oriented source, Q c h = 3.0413 for screw-oriented source). For a Frank-Read source of length 
/ parallel to the x axis and centered at ( x s ,y s ), the source region is approximately bounded 
by an ellipse 


^s d H (x,y) 


{x - x s ) 2 (y - y s ) 2 < 1 


(48) 


(ail) 2 ( a 2 l) 2 

where a± and 02 are calculated to be 2.4610 and 2.2488, respectively. 

With t c , t nvLC and Dg d determined by Eqs. (46), (47) and (48), the Frank-Read source 
is incorporated into the single-slip plane model in the following sense. When the resolved 
shear stress r > r c , a Frank-Read source keeps changing the value of </> 2d at a speed of l/t nuc 
such that dislocation loops enclosing area fl 2d are continuously generated into the system. 
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Mathematically, the source term s 2d in the single-slip plane continuum model in Eq. (8) is 
given by 



(49) 


S' 


where H(z) is the Heaviside function that equals 1 when z > 0 and 0 otherwise, and Xqm is 
the characteristic function in fi 2d , i.e. xqm is 1 hi ^s d an d vanishes elsewhere. 

2.8.2. Incorporation of Frank-Read sources into the three-dimensional continuum model 

Now we derive the expression of the Frank-Read sources in the three-dimensional contin¬ 
uum model. We first write down a formulation for individual Frank-Read sources in three 
dimensions by generalizing the single-slip plane model reviewed above, and then derive a 
formulation based on a source continuum for the three-dimensional continuum model. The 
formulation is affiliated with a slip system, and we omit the slip system superscript A for 
simplicity of notations. 

We first write down a formulation for individual Frank-Read sources in three dimensions. 
Following the construction of three-dimensional dislocation continuum described in Sec. 2.2, 
we generalize the source region fi 2d in the single slip plane model given in Eq. (48) to three 
dimensions by a cylinder Hg d = fi 2d x [—d s i/ 2 , d s i/ 2 ]m, i.e. the intersection of the cylinder 
Hg d with a slip plane is H 2d , and the height of the cylinder is d s \ in the slip plane normal 
direction m, where d s \ is recalled to be the averaged slip plane spacing. 

Since the physical dimension of the source region is much smaller compared to that 
of the domain size D, we further approximate the source region Hg d by a regularized point 
source with same volume at the continuum level, i.e. Xn3 d (') ~ l^s d |^reg(')> where <5 reg (-) 
is a regularized Dirac function of the source region and |Dg d | = 7raia 2 / 2 dsi is the volume 
of flg d . When there are S Frank-Read sources operating, a preliminary way to incorporate 
them into the continuum model in Eq. (42) is to express the source term s as the sum of 
contributions from all the individual Frank-Read sources located at r^, k — 1, 2 , ■ ■ ■ , S: 



s 
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However, in the continuum model, we do not want to resolve individual Frank-Read 
sources. For this purpose, we introduce a source continuum to approximate the collective 
effect of all the individual Frank-Read sources. This is achieved by the following source term 
in the continuum model in Eq. (42): 

s = g{r){r - sign(r)r 0 (r))i7(|r| - r 0 (r)), (51) 


where at any point r, To(r) is the source activation stress at r, and (r — To(r))g(r) measures 
the rate of plastic shear slip initiated at r by the source continuum. 

Now we determine the functions 7o(r) and g( r) from the discrete source model in Eq. (50). 
When the shear stress r exceeds the critical stresses of all the sources that influence point 
r, all the Heaviside functions in Eqs. (50) and (51) can be dropped, a comparison between 
them gives 


g( r) = -vrm g aia 2 6 2 d s i ^ -^<5, 


-'regV 


(52) 


and 


. . Ulg-UCil^b^ d s \ 

T|)(r) = -a*) - -? 


CfSre g(r - r* 


Q 


k 

ch 


log ( — ) ) , where g( r) ^ 0. (53) 


Here we have used the expression of r c in Eq. (46). The source continuum formula of s in 
Eq. (51) with these expressions of To(r) and g{ r) still provides a good approximation to Si n d 
in Eq. (50) when the shear stress r does not exceed the critical stresses of all the sources 
that influence point r, and it becomes exact again when the shear stress r falls below all the 
critical stresses of the sources that influence point r (In this case, s(r) = s in d(r) = 0). 

The functions To(r) and g(r) for the source continuum can be calculated from the ar¬ 
rangement of the discrete Frank-Read sources. In this paper, the locations and parameters 
of the Frank-Read sources are given initially and do not change in the simulation, and we 
calculate these two functions only once for the initial distribution of the sources. 

An example of the Frank-Read source continuum is given in Fig. 3. Fig. 3(a) shows a set of 
individual Frank-Read sources with their positions and lengths generated randomly following 
a normal and a uniform distribution, respectively, and the profiles of its corresponding g( r) 

on several selected slip planes are drawn in Fig. 3(b). It can be observed that g( r) attains 
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Figure 3: Converting a number of Frank-Read sources to a source continuum in a cuboid sample, (a) The 
discrete sources, (b) The corresponding function g( r) in the source continuum obtained by Eq. (52) on 
several slip planes. The unit of the color bar is m g b. 


a relatively high value around the cuboid center, because the number density of individual 
sources are high in the same place in Fig. 3(a). 


2.9. Summary of the governing equations in our continuum model 

To summarize, the derived continuum model based on the DDPFs is constituted mainly 
by the following equations. 


1 . A constitutive stress rule: Given a dislocation substructure described by <ft x and f> x , 
the long-range stress field cr is determined by solving 

= 2/i fsym(Vu) + + £ ^sym(b A ® Vt/; A ) j (54) 

and the equilibrium condition 


V ■ cr = 0, 


(55) 


with boundary conditions given in Eqs. (33) and (34). 

2. A plastic flow rule: The motion of dislocations belonging to the A-th slip system is 
described by an evolution equation for <f x 


d<p x | A | V0 A x X7f> 
~dT + Vn |V^ A | 


AI 


= g\ r) ( t a 


sign(r A )T A (r))tf(T-T U (r)), 


(56) 
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where the dislocation speed v A is determined by the mobility law 


V n = ™g & Votal = m gh X ("hong + T self) > ( 57 ) 

and the shear stress component resolved in the A-th slip system is calculated 
from the long-range stress held cr by Eq. (38), r A lf is the contribution due to the local 
dislocation line tension effect given in Eq. (41), and g x ( r) and r^(r) are two functions 
for the Frank-Read source continuum given by Eqs. (52) and (53), respectively. 

In addition, the Nye dislocation density tensor a, the total dislocation density in the 
system p tot , and the total plastic strain rate ef ot are expressed in terms of the DDPFs 0 and 
0 by Eqs. (22), (44), and (45), respectively. 


2.10. Free energy 

In our DDPF-based continuum model, we can write a total free energy of the system as 

£ ^elastic T ^self; (58) 


where Elastic is the elastic energy 


£ 


elastic 


/.H? 


kXuX 


(b x y 




c-.vu + J2 |^sym(b A ® V0 A ) 


dV— 


/ t-udS 1 

'n t 

(59) 


and £ se if is the energy due to the dislocation line tension effect 


f seif= y ] [ 

, JVL 


p|V'0 A x V0 > 
An 


1 + 


b A • V<^ 


log- 


1 — v (6 A |m A x V0 A |) 2 27rr c |m x V0 


dV. (60) 


This total free energy in Eq. (58) in its form agrees with existing theories using other 
representations of dislocation densities (e.g. Nelson and Toner, 1981; Berdichevsky, 2006; 
Le and Guenther, 2014) and our continuum dislocation dynamics model in a single slip 
plane (Xiang, 2009). 

Since the DDPFs 0 A for all A stay constant during the deformation process, we can write 
£ — £ (u, 0 1 , 0 2 , • • • ,0 A ,---). It can be calculated that inside D, we have 


6 £_ 

hu 


— —V • cr, 
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and 


IW*I i A \ , fi9 , 

£0A ^A (dong + T self) ) (62) 

Recall that the resolved shear stress Tj A is from the long-range stress field cr given by 
Eq. (38), and r s A lf is the contribution due to the local dislocation line tension effect given in 
Eq. (41). Note that the variation of £ se if with respect to 0 A gives the correct leading order 
contribution in r s A lf . 

Using these variations of the free energy, the governing equations in our continuum model 
can be written as 


0 



dt 


6£_ 

5u 


, 8S 

5(j) x 


+ s' 


(63) 

(64) 


The first equation gives the equilibrium condition in Eq. (27), which together with the 
constitutive relation in Eq. (31) and the boundary conditions in Eqs. (33) and (34) form the 
elasticity system. (These boundary conditions can also be obtained from the variation of the 
elastic energy in Eq. (59).) The second equation describes the plastic flow rule (dynamics 
of dislocations), where f^x = —jtx can be understood as the configurational force and 
L x = "' g(6 is the kinetic coefficient. These equations are consistent with the 

governing equations in the DDD models (Hirth and Lothe (1982) and those references in 
the introduction). 

When the system evolves to its equilibrium state, SS/5(p x = 0, which gives rise to the 
micro force balance state Tj A ng + r s A lf = 0. This together with Eq. (63) agree with the models 
for equilibrium states of distributions of straight dislocations (e.g. Berdichevsky, 2006; Le 
and Guenther, 2014). 


3. Numerical implementation of the continuum model 

In this section, we briefly discuss the numerical implementation of our continuum model 
as summarized in Sec. 2.9. We will focus on solving for the long-range stress field cr from 
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Eqs. (54) and (55) with boundary conditions and the plastic flow rule described by the 
evolution of DDPF 0 A in Eq. (56). 

In the simulations in this paper, the computational domain is chosen to be a cuboid 
= [—D/2,D/2] x [—D/2,D/2] x [— L/2,L/2]. The boundary conditions are imposed as 
follows. On the bottom surface, no displacement is allowed along the z direction, which 
is the loading direction, and on the top surface u^\ z=L , 2 = Ug(t) is imposed as a result of 
compression. The shear force is set to be free on these two surfaces. On the other four side 
surfaces, traction free boundary conditions are imposed. The DDPF in this paper takes 
the form in Eq. (13) with inter-slip plane distance d s \ = 1005, and does not change in the 
plastic deformation process. 

3.1. Finite element formulation for solving for the long-range stress field 

The long-range stress field satisfying Eqs. (54) and (55) subject to boundary conditions 
in Eqs. (33) and (34) yields a weak form that f n a i ong : (Vv)dE = f m t b • vd S for any test 
vector function v E {v|v = 0, on <9h2d}- Replacing the stress field by the constitutive stress 
rule given by Eq. (54), we obtain the weak form for the displacement field u 

f Vv : C : VudE = / t • vd S — 2ji Y'' f </> A sym(b A ® V'0 A ) : (Vv)dV, (65) 
J J <9f2t ^ J n 

where the forth-order tensor C is defined in Eq. (32). 

In our simulations, is meshed by C3D8 bricks. We then discretize the weak form in 
Eq. (65) to get a linear algebraic equation system as 

ApE u FE = fFE) (66) 

where Ufe is a vector of 31V dimensions containing all nodal values of u with N being 
the total number of nodes, ARe is known as the stiffness matrix, and fpp is assembled by 
discretizing the right hand side of Eq. (65). 

If the last term is omitted, Eq. (65) is the weak formulation that is widely used in the FE 
formulation in linear elasticity. Hence many tools well-developed for solving purely elastic 

problems, such as meshing, assembling and inversion of A'fe, can be inherited by the FE 
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formulation proposed here. The contribution from 0 A and in Eq. (66) can be treated the 
same way as a body force to the linear system. 


3.2. Finite difference scheme for the evolution of DDPF 0 A 

We discretize the evolution equation of DDPF <f x in Eq. (56) using a finite difference 
scheme. Here the grid points of <f x are chosen coinciding with the vertices of the C3D8 
bricks. Combined with the mobility law in Eq. (57), Eq. (56) can be written as 

d(j) x 


dt 


+ ( T tone + T ilf) l m> x V| A X | = s> 


(67) 


The temporal derivative of <f x is approximated by the forward Euler scheme. For the spatial 
derivatives of </> A , we follow the methods in Xiang et al. (2003) that the first-order upwind 
scheme is used for those in the term associated with the long-range stress field r longl mA X 
V0 A |, and the central difference scheme is used to calculate those spatial derivatives in the 


term r A lf | 


m 


x V0 A | where r s A lf is given by Eq. (41). 


Moreover, the regularized ^-function in the source term s A (Eqs. (51)-(53)) is given by 


^reg(r) — 


7r 


Asf 7T 2 — 4 


7r|m A x r 


cos 


Asi 


+ 1 


2A So 


7T( m • r 

cos- 7 -- + 1 


As 2 


( 68 ) 


for all r E {r||m A x r| < Asi, |m A • r| < AS 2 }, where Asi and As 2 are two smoothing 


parameters. 


4. Numerical examples 

In this section, the derived continuum model is validated through comparisons with DDD 
simulations. All results presented are obtained by using 10 x 10 x 20 C3D8 bricks in the 
finite element discretization of the cuboid simulation domain described at the beginning of 
Sec. 3. 

4-1. A single Frank-Read source under a constant applied strain 

This example is aimed to provide an illustration of the continuum model. For the cuboid 
simulation cell L = 240006 and D = 96006. A Frank-Read source of length / = 4006 








with activation stress 7.8 x 10~ 4 /r is placed at the center of 0, and a 0.3% constant strain 
is applied by compression on the top surface. All simulations start with a dislocation-free 
state. 



Figure 4: Snap shots of the distribution of dislocation curves generated by an operating Frank-Read source 
that locates at the center of the simulation domain ft. The source releases dislocation loops in response to 
a constant strain by compression on the top surface of ft. The length unit of ft is its height L. 

In Fig. 4, we plot the contour lines of the DDPF (j) on one of the slip planes, which give 
rough locations of the dislocation curves. It is observed that in response to the applied strain, 
the source keeps releasing dislocation loops, which exit hi from its side surfaces. As a result, 
the resolved shear stress drops during this loop-releasing process, and so does the pressure 
on the top surface as shown in Fig. 5. These findings agree with the common understanding 
about the role played by a Frank-Read source: It releases dislocation loops so as to soften 
the materials. When the time t is about 14000 L/(m g /j,b), the resolved shear stress finally 
falls below the source activation stress, the source is then deactivated, see Fig. 5. 

The displacement vector u on the domain boundary dfl is shown in Fig. 6. It can be 
seen that the orientations of the displacement gradients are in the normal direction of the 
activated slip plane that contains the Frank-Read source, and they are localized near this 
slip plane. This is because the dislocation loops that have left the specimen form surface 
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Figure 5: As the Frank-Read source keeps releasing dislocation loops, both the resolved shear stress at 
the source and the pressure on the top surface of the simulation domain drop. When the time t is about 
14000L/(m g ^&), the resolved shear stress falls below the source activation stress (indicated by the dashed- 
dotted line), and the source is thus deactivated. 



(a) ui 


(b) u 2 


(c) u 3 


Figure 6: The three components of the displacement vector u on the simulation cell surface dCl. The length 
unit is L , the height of 17. 


steps on the activated slip plane (in a smooth sense). 


4-2. Comparison with DDD simulations 

To further validate the continuum model, we compare its numerical results with the DDD 
simulation results obtained by El-A wady et al. (2008). Following their DDD simulations, 
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we choose a micro-pillar of nickel. The loading axis is < 001 >, and a single slip system 
is activated with slip direction [Oil] and slip normal (111). The Schmid’s factor m s is thus 
0.4050. The shear modulus is 76GPa, the Poisson’s ratio is 0.31, and the length of the 
Burgers vector is b = 0.25nm. The strain rate is 200s” 1 . The dislocation glide mobility 
m g is unspecified by El-Awady et al. (2008), we here follow Senger et al. (2008) to let 
m g = 10 4 /(Pa-s). In our simulations, the micro-pillars are chosen to be cuboids as described 
at the beginning of Sec. 3 for the ease of the finite element implementation. The sample 
sizes here are defined to be the size of the cuboid base D. The (height to base size) aspect 
ratio of the micro-pillars is L/D = 3. We choose D = 0.5/im or lprn. 


Sample 

Mean source 

length (pm) 

Standard 

deviation (pm) 

Max source 

length (pm) 

T min 

(MPa) 

Flow stress 

(MPa) 

1 

0.1835 

0.1104 

0.4832 

170.7 

420.0 

2 

0.2403 

0.0935 

0.4423 

132.8 

331.5 

3 

0.2039 

0.0992 

0.3340 

145.5 

357.9 

4 

0.1928 

0.1083 

0.3711 

158.9 

388.8 


Table 1: Statistics of the individual sources whose lengths are obtained randomly following a uniform 
distribution within [20nm, D] for samples of size D = 0.5pm. The locations of the sources are also determined 
randomly following the uniform distribution over the micro-pillar. 


Sample 

Mean source 

length (pm) 

Standard 

deviation (pm) 

Max source 

length (pm) 

T° ■ 

1 nun 

(MPa) 

Flow stress 

(MPa) 

1 

0.4066 

0.2252 

0.9196 

83.8 

196.4 

2 

0.4047 

0.2412 

0.9995 

91.0 

230.2 

3 

0.4596 

0.2333 

0.9605 

80.7 

204.0 

4 

0.4709 

0.2275 

0.9349 

80.3 

201.2 


Table 2: Same as Table 1 for samples of size D = lpm. 


In our simulations, the lengths of the Frank-Read sources are generated randomly fol¬ 
lowing a uniform distribution within [20nm, D\. The initial dislocation density po is also 
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randomly generated within the range 1.6 ~ 4 x 10 12 m -2 , and these pre-existing dislocation 
segments of the sources are uniformly assigned to the twelve slip systems of fee nickel. The 
statistics of the initial source distributions for various samples are listed in Tables 1 and 2. 
With these isolated Frank-Read sources generated, the corresponding source continuum is 
developed following Eqs. (51)-(53). For the source character parameter C s in the critical 
stress formula in Eq. (46), we choose it to be 1.35 and 2.02 for a single edge source and a 
screw source, respectively, to agree with the simulation set-up in El-Awady et al. (2008). 

In Tables 1 and 2, we also show the values of r)h n , which is defined as the minimum 
value of the activation stress T 0 (r) of the Frank-Read source continuum given by Eq. (53), 
over the points where g( r) ^ 0 inside the pillar. The minimum activation stress r^ in will be 
used in the derivation of the scaling law in the next section. 
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Figure 7: Stress-strain curves obtained by the simulations of our continuum model. The vertical bars denote 
the ranges of the flow stress predicted by the DDD simulations of El-Awady et al. (2008), whose averaged 
values are shown by the black dots. 


The stress-strain curves obtained by our continuum model are shown in Fig. 7 for samples 
of sizes D = 0.5/nn and D = 1/irn. Good agreement with the DDD results by El-Awady 
et al. (2008) are observed. Firstly, both methods predict an initially elastic regime and an 
almost perfectly plastic regime, where work-hardening effect is barely observed. Secondly, 
both methods suggest that the engineering stress stays roughly unchanged or oscillate around 
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some constant value in the regime of perfect plastic deformation, and this stress is measured 
as the flow stress of the micro-pillars. The “smaller-being-stronger” size effect on crystalline 
strength, which is indicated by the flow stress, is observed in both our simulations and 
in the DDD simulations of El-Awady et ah (2008). Moreover, statistical effects in the 
flow stress are seen in the simulations of both methods. Such statistical effects have also 
been examined in other literature (e.g. El-Awady et al., 2009; hi et ah, 2014). To make 
quantitative comparisons between the results of our model and those of the DDD simulations 
of El-Awady et ah (2008), the respective ranges and averaged values for the flow stress 
recorded by El-Awady et ah (2008) are drawn in Fig. 7. It can be seen that the flow stresses 
calculated by our continuum model all fall into the respective ranges predicted by their DDD 
simulations. 

These comparisons show that our continuum model is able to provide an excellent sum¬ 
mary of the corresponding underlying dynamics of discrete dislocations. In the next section, 
we will use it to quantitatively study the size effect on strength observed in the uniaxial 
compression tests of micro-pillars. 

5. Size effect on strength of single-crystalline micro-pillars 

5.1. Comparison with the experimental data 

Now we investigate the “smaller-being-stronger” size effect in micro-pillars using our 
continuum model. For the initialization of the Frank-Read sources, we follow Shishvan and 
Van der Giessen (2010) that, in analogy to the distribution of grain sizes in polycrystals 
which has been experimentally measured, the Frank-Read source size l follows a log-normal 
distribution with the probability density function 

1 (logi-logi m ) 2 

P (0 = —f= —7 e ^ (69) 

v 2vrcr s d< 

with two parameters l m and <r s d to be determined. The parameter l m which can be considered 
as the effective mean source length should decrease with the pillar size D. The locations 
of these source segments follow the uniform random distribution over the micro-pillar. The 
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single-ended sources (Parthasarathy et al., 2007) are included, and this happens when part 
of the source segments are outside the pillar. 

To validate the initial source distribution in Eq. (69) and determine the parameters, the 
numerical results by using our continuum model are compared with the experimental data 
of nickel micro-pillars by Dimiduk et al. (2005). In our simulations, most parameters are 
chosen the same as those used in Sec. 4.2 with the following exceptions in accordance with 
Dimiduk et al. (2005): the loading axis is set along [269] and the singly active slip system is 
of the slip direction [101] and slip normal (111), the Schmid factor is m s = 0.48, the shear 
modulus is 78GPa, and the aspect ratio of the pillar ( L/D ) is chosen randomly between 2 
and 3. The total density of all source segments is 3 x 10 -12 m -2 following Dimiduk et al. 


(2005). 


To fit their experimental results, we choose the effective mean source length / m = a m D , 
where a m = 1/15. The standard deviation a s d is determined under the assumption that the 
probability of a source segment that is longer than D is no more than 10~‘, which gives 
<j s d = 0.4. The computed stress-strain curves with these values of parameters are shown in 
Fig. 8 for different samples varying in size. It can be seen that the obtained flow stresses 
are in good agreement with the experimental data of Dimiduk et al. (2005). The size effect 
on strength is also clearly observed from our simulation results. 

5.2. Scaling law of the size effect on micro-pillar strength 

To describe how the pillar strength depends on its size D, we propose a formula 



(70) 


where cr[[ ow is some constant stress, m s is recalled to be the Schmid factor, and a e Q is a 
dimensionless constant such that a c ^D measures the length of the weakest source in pillar 
of size D. Eq. (70) suggests a scaling law 



(71) 


A comparison of Eq. (70) with experimental data (to be shown at the end of this subsection) 

suggests that Eq. (70) pocesses a wider effective zone, that is, it is valid for many (at least 
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Figure 8: Stress-strain curves obtained by simulations using our continuum model with initial arrangement 
of Frank-Read sources whose length distribution following the log-normal one given by Eq. (69). Four 
simulations are performed for each size of the micro-pillar D that ranges from 1/rm to 20/rm. The vertical 
bars identify the ranges of the experimentally measured flow stress by Dimiduk et al. (2005). 


35 









































three) types of f.c.c. pillars of size ranging from submicrons to tens of microns. We now 
rationalize Eq. (70) by taking the following two steps. 
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Figure 9: (a) Linear dependence between the flow stress (Jflow and the minimum activation stress of the 
source continuum r^ in given in Eq. (72) (the dashed line), with fitted parameter er){ ow . (b) The dependence 
of the minimum activation stress rT n on the sample size D given in Eq. (74) (the dashed line), where 
= 1/6 and r c = 0.66. These relations are validated by numerical results obtained using the continuum 
model for pillars with different sizes (specified in the legend) shown by the dots in (a) and (b). 


We first show that the flow stress is related to the minimum activation stress r„ lin by 

^0 


T ■ 

_ _ min | 0 

^flow i 


m* 


flow’ 


(72) 


where a{] ow is some constant stress, and recall that m s is the Schmid factor. This means 
that the flow stress is determined by the weakest Frank-Read source inside the pillar. Com¬ 
parisons of the results of this linear relation and the computed flow stresses by using the 
continuum model for pillars with different sizes are shown in Fig. 9(a), in which the pa¬ 
rameter Ofl ow is fitted from the numerical results. The excellent agreement seen in these 
comparisons demonstrates that Eq. (72) provides a good quantitative description for the 
flow stress <7fl ow in terms of r° iT1 . 

The next step is to relate r° lin to the sample size D. We follow the formula of the 
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activation stress of a single Frank-Read source to assume 


T ■ 
mm 



( 73 ) 


where C s is dependent on the source character chosen to be 1, l e g can be considered as an 
effective source length and it is assumed to be a fraction of D. Using simulation results of 
the continuum model for pillars with different sizes, fitting r° lin against l e g with reference to 
Eq. (73) gives that l e g = a e gD , where the parameter a e g =1/6. Physically, this means that 


given a micro-pillar of size D , the longest source length is most likely a e gD. Hence r^ in is 
related to D by 



P4) 


This relation is validated by numerical results obtained using the continuum model for pillars 
with different sizes as shown in Fig. 9(b). 

Combining Eqs. (72) and (74), we obtain the scaling law in Eq. (70). Predictions of 
Eq. (70) are in good agreement with the experimental data as shown in Fig. 10 for single 
crystal nickel, aluminium and copper pillars. In each comparison in the figure, the resolved 
flow stress Tfl ow = m s ag ow is plotted against the pillar size D, with a e g = 1/6, r c = 0.66 and 
°flow being a fitted parameter in the scaling law. 

It is noted that Eq. (74) has also been employed to predict the flow stress in polycrys¬ 
talline thin-hlm structures (von Blanckenhagen et al., 2001; Gruber et al., 2008). In their 
works, a e g is suggested to be 1/4 to 1/3. 

5.3. Other behaviors of the micro-pillars 

In this subsection, we discuss other plastic behaviors of micro-pillars that can be captured 
by the continuum model. 

Firstly, recall that in our continuum model, the dislocation networks within the specimen 
(after local homogenization) are described by the DDPFs as defined in Sec. 2.2 and illustrated 
by Fig. 2, that is, the dislocation distribution on each slip plane is described by the contour 
curves of the DDPF 0, and the distribution of the slip plane is represented by another DDPF 
fj. Fig. 11 shows some snapshots of the dislocation substructures in some simulations using 
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(a) Nickel 


(b) Aluminium 


CD 

SI 

CD 



(c) Copper < 111 > 


(d) Copper < 123 > 


Figure 10: The log (D)/D scaling law for the size-dependent effect of the flow stress of a micro-pillar given in 
Eq. (70) resolved in the slip planes is examined by the experimental data for single crystal nickel (Dimiduk 
et ah, 2005), aluminium and copper (Uchic et ah, 2009) micro-pillars. 

the continuum model for pillars of size D = 2.4/mr, 5pm, and 10pm. The distributions of 
dislocation curves in the figure look smoothly varying. This is because the continuum model 
only resolves the dislocation microstructures in an average sense (see Sec. 2.2). 

We can also keep track of the two quantities that are commonly used to describe the 
state of the specimen in a plastic deformation process: the total dislocation density p tot 
given by Eq. (44) and the plastic strain rate ef ot given by Eq. (45). Numerical results of 
Ptot and ej? ot obtained by simulations using our continuum model are plotted in Fig. 12. 
These results along with the stress-strain curves presented in Figs. 7 and 8 suggest that the 
following evolution process may take place inside the micro-pillars when being compressed 
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(a) D = 2.4/Ltm (b) D = 5/mi (c) D = 10^m 

Figure 11: Snapshots of the dislocation substructures (in an average sense as discussed in Sec. 2.2) in some 
simulations using our continuum model for pillars of size (a) D = 2.4/im, (b) 5//m, and (c) 10/rm. The 
snapshots are taken from the perfect plastic deformation regime as shown by the stress-strain curves in 
Figs. 7 and 8. 




(a) (b) 

Figure 12: Evolution of (a) the total dislocation density determined by Eq. (44) and (b) the plastic strain 
rate determined by Eq. (45) in micro-pillars with various sizes, obtained by using the continuum model. 
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under a constant applied strain rate. When the elastic limit of the samples is reached, 
dislocation sources start to release dislocation loops, resulting in plastic flows and a rise in 
the total dislocation density inside the pillars. After a (relatively) short period, the system 
reaches a steady state corresponding to the perfectly plastic regimes in Figs. 7 and 8. At 
this steady state, the applied strain rate is fully accommodated by the dislocation motion 
and the resolved shear stress ceases to increase. Another interesting phenomenon observed 
from Fig. 12(b) is that the values of ej) ot converge to a same value for samples of various size. 
This is because this converged values are determined by the applied strain rates which are 
the same for all samples in the compression tests. 



(a) 2.5% strain (b) 9% strain (c) 12% strain 

Figure 13: Shapes of a deformed micro-pillar of size D = 5/xm under various applied strain: the cuboids 
formed by the dashed-lines describe the original shapes of the pillar. 

By using the continuum model, we are also able to track the shape changes of the 
micro-pillars. Given u the displacement field, r + u is the position of a point, whose initial 
position is at r. Examples of the profile of a deformed pillar during compression are shown in 
Fig. 13. Note that in our simulations, spatial variation in the plastic shear slips induced by 
a distribution of Frank-Read sources has been locally averaged over many discrete sources in 
the source continuum formulation. The pillar profiles in Fig. 13 agree in the averaged sense 

with the experimental observations (e.g. Fig. 4 in Dimiduk et al. (2005)). With smaller 
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number of sources, relatively strongly non-uniform displacement fields along z direction can 
be observed in our simulations, see Fig. 14 for an example, which more accurately agree 
with the slip-band structures observed in the experiments. For an extreme case, if we check 
the surface of a micro-pillar containing only one operating Frank-Read source as shown in 
Fig. 6, localized displacement gradient is clearly observed around the activated slip plane. 
These agree with the conclusion of Akarapu et al. (2010) drawn from simulations using a 
hybrid elasto-viscoplastic model which couples DDD that such localized deformation, which 
is strong for small-size pillars, is due to the heterogeneous dislocation distributions that lead 
to clusters of sources. More efficient numerical implementation method of our continuum 
model with fine meshes will help resolve more accurately the pillar profile observed in the 
experiments, which will be developed in the future work. 





(a) ui (b) u 2 (c) u 3 

Figure 14: Displacement field u = (ui,u 2 ,u 3 ) along a vertical edge of a micro-pillar with size D = lym and 
8 Frank-Read sources under 0.3% applied strain. The length unit is L. 


6. Conclusion 

In this paper, we have presented a dislocation-density-based continuum model to study 
the plastic behavior of crystals, in which the dislocation substructures are represented by 
pairs of DDPFs. A discrete dislocation network in three dimensions is approximated by 
a dislocation continuum after local homogenization of dislocation ensembles within some 
representative volume. For the A-th slip plane in the dislocation continuum, a DDPF o x is 
defined such that </> A restricted on each slip plane describes the plastic slip across the slip 
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plane and identifies the distribution of dislocation curves by its contour lines, and another 
DDPF ip x is employed to describe the slip plane distribution by its contour surfaces. The 
major advantage of this three-dimensional continuum model lies in its simple representa¬ 
tion of distributions of curved dislocations. The connectivity condition of dislocations is 
automatically satisfied with this representation of dislocations. 

Based on the DDPFs, the plastic deformation process of crystals can be formulated 
by a system of equations as given by Eqs. (54)-(57). We have shown that the equation 
system provides an effective summary over the underlying discrete dislocation dynamics, 
including the long-range elastic interaction of dislocations, dislocation line tension effect, and 
dislocation multiplication by Frank-Read sources. Numerically, a finite element formulation 
is proposed to compute the long-range stress held. 

The continuum model is validated by comparisons with DDD simulations and experi¬ 
mental data. As one application of the continuum model characterized by DDPFs, the size 
effect on the strength of micro-pillars is studied, and the pillar how stress is found scaling 
with its (non-dimensionalized) pillar size D by log (D)/D. 

Future work may include generalizations of the continuum model characterized by DDPFs 
to incorporate the out-of-slip-plane dislocation motion (cross-slip or climb) in which the con¬ 
tour surfaces of the DDPF may be curved and evolved in the plastic deformation. We 
will also take into account in the continuum model other short-range dislocation interactions 
such as dislocation reactions and junction formation, as well as dislocation interactions with 
other defects (e.g. Xiang and Srolovitz, 2006; Chen et ah, 2010; Zhu and Xiang, 2014). 
Efficient numerical implementation method will also be considered in the future work. 
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